---
title: "Replication Code for: 'Domestic Costs of Nuclear Deterrence: Voter Turnout and Nuclear Weapons Testing'"
author: "Unislawa M. Williams, Mya Whiles, Tinaz Pavri"
output:
  word_document: default
  pdf_document: default
  number_sections: true
date: "28 April 2025"

---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, results = "markup")
#setwd("")
```



```{r, echo=FALSE, results='asis'}
library(dplyr)
library(fixest)
library(stringr)
library(knitr)
library(kableExtra)
library(car)
library(ggplot2)
library(modelsummary)


# Read in the data
load("discop.Rda") #External data to be downloaded from Marsh (2022). "Replication Data for: 'Trauma and Turnout: The Political Consequences of Traumatic Events'" Harvard Dataverse. https://doi.org/doi:10.7910/DVN/MXKMSE.
compensated <- read.csv("Compensated.csv")
merged_data <- discop %>%
  left_join(compensated, by = c("state", "county", "year"))

```


```{r, echo=FALSE, results='asis'}
#This code reproduces Table 1

years <- c(1980, 1984, 1988, 1992, 1996)

merged <- merged_data %>%
  filter(
    state %in% c("UT", "AZ", "NV"),
    year %in% years,
    !is.na(`Compensated1990`)
  ) %>%
  mutate(
    treated = Considered_For_Compensation_But_Not_Compensated1990,
    ratio_gop_dem = votes_gop / votes_dem,
    perc_bipoc = (pop_tot - pop_wh) / pop_tot,
    pop_tot = pop_tot,
    post = post
  )

modelUT_1988_1992 <- feols(
  turnout ~ treated * post + pop_tot + income + perc_bipoc + ratio_gop_dem | county + year,
  data = merged %>% filter(year %in% c(1988, 1992), state %in% "UT")
)

#this analysis uses an older version of fixest package

modelUT_1988_1996 <- feols(
  turnout ~ treated * post + pop_tot + income + perc_bipoc + ratio_gop_dem | county + year,
  data = merged %>% filter(year %in% c(1988, 1996), state %in% "UT")
)

model_1988_1992 <- feols(
  turnout ~ treated * post + pop_tot + income + perc_bipoc + ratio_gop_dem | county + year,
  data = merged %>% filter(year %in% c(1988, 1992))
)

model_1988_1996 <- feols(
  turnout ~ treated * post + pop_tot + income + perc_bipoc + ratio_gop_dem | county + year,
  data = merged %>% filter(year %in% c(1988, 1996))
)

models <- list(
  "1988 vs. 1992: UT" = modelUT_1988_1992,
  "1988 vs. 1996: UT" = modelUT_1988_1996,
  "1988 vs. 1992: AZ/NV/UT" = model_1988_1992,
  "1988 vs. 1996: AZ/NV/UT" = model_1988_1996
)

modelsummary(
  models,
  output = "markdown", 
  statistic = "({std.error})", 
  coef_map = c(
    "treated:post" = "RECA (1990): Treated x Post", 
    "pop_tot" = "Total Population", 
    "income" = "Income", 
    "perc_bipoc" = "BIPOC Voters (%)", 
    "ratio_gop_dem" = "GOP:Dem Ratio"
  ),
  stars = c('.' = 0.1, '*' = 0.05), 
  gof_omit = "AIC|BIC|Log.Lik"
)

# Create a word table
modelsummary(
  models,
  output = "markdown",
  statistic = "({std.error}) {p.value}", 
  stars = c('.' = 0.1, '*' = 0.05, '***' = 0.01), 
  gof_omit = "AIC|BIC|Log.Lik" 
)
```

```{r, echo=FALSE, results='asis'}
#This code supports hypothesis 2
######How many voters?###############################################################
aggregate_turnout_UT <- merged %>%
  filter(
    year %in% c(1992), 
    state == "UT",
    Considered_For_Compensation_But_Not_Compensated1990 == 1) 
voter_increase_UT = .104 * sum(aggregate_turnout_UT$pop_tot) #.1 is the coeff on DID
print(voter_increase_UT)
#Turnout increase in UT in 1992
#and
#Turnout increases in UT/NV/AZ in 1992 and 1996
aggregate_turnout <- merged %>%
  filter(
    year %in% c(1992), 
    Considered_For_Compensation_But_Not_Compensated1990 == 1) 
voter_increase_92 = .091 * sum(aggregate_turnout$pop_tot) #.091 is the coeff on DID
print(voter_increase_92)

aggregate_turnout <- merged %>%
  filter(
    year %in% c(1996), 
    Considered_For_Compensation_But_Not_Compensated1990 == 1) 
voter_increase_96 = .085 * sum(aggregate_turnout$pop_tot) #.085 is the coeff on DID
print(voter_increase_96)
print(c("Total additional votes in the 1992 and 1996 elections:", sum(voter_increase_92, voter_increase_96))) 
#This number is likely an overestimate as total population is not the same as the population of elgible voters, i.e. children etc.  However, it highlights the aggregate size of the effect.
#According to Dept of Justice Downwinder total pending, approved and denied claims in 2022 included for all affected areas was
#Total 30,066 RECA Downwinder Claims (Congressonal Report 2022)
```

The total of additional votes in the 1992 and 1996 elections across the the three states was likely lower than the estimate of 94,974 voters.  This is likely an overestimate as total population is not the same as the population of eligible voters, i.e. children etc.  However, it highlights the aggregate size of the effect. According to Dept of Justice, Downwinder total pending, approved and denied claims in 2022 for all affected areas were lower than the 1992 increase. Total 30,066 RECA Downwinder Claims (Congressional Report 2022) were reported.  



```{r, echo=FALSE, results='asis'}
######Testing Assumptions############################################################
#1. Plotting the trends shows some diff but not large
pre_trends <- merged %>%
  filter(year < 1990, state == "UT") %>%
  group_by(year, treated) %>%
  summarise(mean_turnout = mean(turnout, na.rm = TRUE), .groups = "drop")


ggplot(pre_trends, aes(x = year, y = mean_turnout, color = factor(treated), group = factor(treated))) +
  geom_line(size = 1) +
  geom_point(size = 3) +
  labs(
    title = "Parallel Trends Assessment: Turnout Pre-Treatment (UT: 1980-1988)",
    x = "Year",
    y = "Mean Turnout",
    color = "Treated Group"
  ) +
  theme_minimal()

pre_trends <- merged %>%
  filter(year < 1990) %>%
  group_by(year, treated) %>%
  summarise(mean_turnout = mean(turnout, na.rm = TRUE), .groups = "drop")


ggplot(pre_trends, aes(x = year, y = mean_turnout, color = factor(treated), group = factor(treated))) +
  geom_line(size = 1) +
  geom_point(size = 3) +
  labs(
    title = "Parallel Trends Assessment: Turnout Pre-Treatment (UT: 1980-1988)",
    x = "Year",
    y = "Mean Turnout",
    color = "Treated Group"
  ) +
  theme_minimal()

#2. The placebo model testing whether 1988 or 1984 in UT or all three states were insignificant
# 1984 vs. 1988
placebo_data <- merged %>%
  filter(year %in% c(1984, 1988)) %>%
  mutate(post = ifelse(year == 1988, 1, 0))  # Define post for placebo as 1988

# Run placebo DiD model
placebo_model <- feols(
  turnout ~ treated * post + pop_tot + income + perc_bipoc + ratio_gop_dem | county + year,
  data = placebo_data %>% filter(year %in% c(1984, 1988))
)
placebo_model_UT <- feols(
  turnout ~ treated * post + pop_tot + income + perc_bipoc + ratio_gop_dem | county + year,
  data = placebo_data %>% filter(year %in% c(1984, 1988), state %in% "UT")
)
# There is very little data for 1980 
placebo_data <- merged %>%
  filter(year %in% c(1980, 1984)) %>%
  mutate(post = ifelse(year == 1984, 1, 0))  # Define post for placebo as 1988

# Run placebo DiD model
placebo_model1984 <- feols(
  turnout ~ treated * post + pop_tot + income + perc_bipoc + ratio_gop_dem | county + year,
  data = placebo_data %>% filter(year %in% c(1980, 1984))
)
placebo_model_UT1984 <- feols(
  turnout ~ treated * post + pop_tot + income + perc_bipoc + ratio_gop_dem | county + year,
  data = placebo_data %>% filter(year %in% c(1980, 1984), state %in% "UT")
)

placebo_models = list(
  placebo_model,
  placebo_model_UT,
  placebo_model1984,
  placebo_model_UT1984)

modelsummary(
  placebo_models,
  output = "markdown", 
  statistic = "({std.error})", 
  coef_map = c(
    "treated:post" = "RECA (1990): Treated x Post", 
    "pop_tot" = "Total Population", 
    "income" = "Income", 
    "perc_bipoc" = "BIPOC Voters (%)", 
    "ratio_gop_dem" = "GOP:Dem Ratio"
  ),
  stars = c('.' = 0.1, '*' = 0.05), 
  gof_omit = "AIC|BIC|Log.Lik"
)
```
The above models show placebo effects for 1984 and 1988 do not have statistical significance across a range of specifications.

```{r, echo=FALSE, results='asis'}
#This code assesses the results in Table 1 based on different categorizations of counties

years <- c(1980, 1984, 1988, 1992, 1996)

merged <- merged_data %>%
  filter(
    state %in% c("UT", "AZ", "NV"),
    year %in% years,
    !is.na(`Compensated1990`)
  ) %>%
  mutate(
    treated = NotCompensated_with_Clark, 
    ratio_gop_dem = votes_gop / votes_dem,
    perc_bipoc = (pop_tot - pop_wh) / pop_tot,
    pop_tot = pop_tot,
    post = post
  )
#for different specification of treated counties use 
#NotCompensated_with_Washington
#NotCompensated_with_Mohave
#NotCompensated_without_Coconino
#NotCompensated_with_Clark


modelUT_1988_1992 <- feols(
  turnout ~ treated * post + pop_tot + income + perc_bipoc + ratio_gop_dem | county + year,
  data = merged %>% filter(year %in% c(1988, 1992), state %in% "UT")
)

#this analysis uses an older version of fixest package

modelUT_1988_1996 <- feols(
  turnout ~ treated * post + pop_tot + income + perc_bipoc + ratio_gop_dem | county + year,
  data = merged %>% filter(year %in% c(1988, 1996), state %in% "UT")
)

model_1988_1992 <- feols(
  turnout ~ treated * post + pop_tot + income + perc_bipoc + ratio_gop_dem | county + year,
  data = merged %>% filter(year %in% c(1988, 1992))
)

model_1988_1996 <- feols(
  turnout ~ treated * post + pop_tot + income + perc_bipoc + ratio_gop_dem | county + year,
  data = merged %>% filter(year %in% c(1988, 1996))
)

models <- list(
  "1988 vs. 1992: UT" = modelUT_1988_1992,
  "1988 vs. 1996: UT" = modelUT_1988_1996,
  "1988 vs. 1992: AZ/NV/UT" = model_1988_1992,
  "1988 vs. 1996: AZ/NV/UT" = model_1988_1996
)

modelsummary(
  models,
  output = "markdown", 
  statistic = "({std.error})", 
  coef_map = c(
    "treated:post" = "RECA (1990): Treated x Post", 
    "pop_tot" = "Total Population", 
    "income" = "Income", 
    "perc_bipoc" = "BIPOC Voters (%)", 
    "ratio_gop_dem" = "GOP:Dem Ratio"
  ),
  stars = c('.' = 0.1, '*' = 0.05), 
  gof_omit = "AIC|BIC|Log.Lik"
)


```

The above models show results are robust to interpretations of county classifications. 